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^-T" ' Abstract 

'NT ■ 

^— N . We analyze the periodically forced Kuramoto model. This system consists 

f^ , of an infinite population of phase oscillators with random intrinsic frequencies, 

global sinusoidal coupling, and external sinusoidal forcing. It represents an 
$H ' idealization of many phenomena in physics, chemistry and biology in which 

mutual synchronization competes with forced synchronization. In other words, 
the oscillators in the population try to synchronize with one another while also 
trying to lock onto an external drive. Previous work on the forced Kuramoto 
model uncovered two main types of attractors, called forced entrainment and 
mutual entrainment, but the details of the bifurcations between them were un- 
clear. Here we present a complete bifurcation analysis of the model for a special 
case in which the infinite-dimensional dynamics collapse to a two-dimensional 



system. Exact results are obtained for the locations of Hopf, saddle-node, and 
Takens-Bogdanov bifurcations. The resulting stability diagram bears a striking 
resemblance to that for the weakly nonlinear forced van der Pol oscillator. 

Abbreviated title: Forced Kuramoto model 



1 INTRODUCTION 

The study of synchronization is a classic topic in nonUnear science. 
Sometimes the concern is with mutual synchronization, as in Huy- 
gens's 1665 discovery of the sympathy of pendulum clocks. In other 
situations, one is more interested in forced synchronization, as in the 
injection locking of a laser or the entrainment of circadian rhythms by 
the daily light-dark cycle. Here w^e consider a simple mathematical 
model in which both types of synchronization are present simulta- 
neously, creating a conflict between them. What happens w^hen a 
network of dissimilar but mutually coupled oscillators is also driven 
by an external periodic force? For a natural generalization of the Ku- 
ramoto model, the interaction of forcing, coupling, and randomness 
leads to a rich set of collective states and bifurcations. We explain all 
of these phenomena analytically, using an ansatz recently introduced 
by Ott and Antonsen. 

1 Introduction 

In 1975 Kuramoto proposed an elegant model for an enormous population of 
coupled biological oscillators [Kuramoto 1975, 1984]. Each oscillator was de- 
scribed solely by its phase, with amplitude variations neglected; the oscilla- 
tors were coupled all-to-all, with equal strength; the interaction between them 
was purely sinusoidal, with no higher harmonics; and their intrinsic frequen- 
cies were randomly distributed across the population according to a symmet- 
ric bell-shaped distribution. All of these simplifying assumptions helped Ku- 
ramoto make headway on what would otherwise have been a hopelessly in- 
tractable many-body, nonlinear dynamical system. By means of an ingenious 
self-consistency argument, he was able to show analytically that the system 
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could undergo a phase transition to mutual synchronization, once the coupling 
between the oscillators exceeded a certain threshold. 

Over the past three decades, many researchers have shed light on the math- 
ematical aspects of collective synchronization by studying Kuramoto's model 
and its close relatives [Strogatz 2000, Acebron et al. 2005]. And, somewhat sur- 
prisingly in view of its simplicity, the model has also been shown to be relevant 
to a variety of physical systems [Pikovsky et al. 2001, Strogatz 2003]. Examples 
range from electrochemical oscillators [Kiss et al. 2002, 2008] and Josephson 
junction arrays [Wiesenfeld et al. 1996] to coupled metronomes [Pantaleone 
2002], collective atomic recoil lasing [von Cube et al. 2004], and neutrino flavor 
oscillations [Pantaleone 1998]. 

One way to extend the model is to allow for the effects of external forcing. 
This generalization is theoretically natural, but it is also motivated in part by 
experimentally observed phenomena [Kiss et al. 2008]. For example, consider 
the way that the daily cycle of light and darkness helps to entrain our sleep, 
body temperature, and other circadian rhythms to the world around us [Moore- 
Ede et al. 1982, Dunlap et al. 2003, Foster and Kreitzman 2005]. Like all 
mammals, each of us has a circadian pacemaker, a network of thousands of 
specialized clock cells located in the region of the hypothalamus known as the 
suprachiasmatic nuclei, just above where the optic nerves criss-cross as they 
make their way back to the brain. These cells have been shown experimentally 
to be intrinsically oscillatory [Welsh et al. 1995] and their distribution of natural 
frequencies has been measured [Liu et al. 1997]. The pacemaker cells are 
also known to be mutually coupled, though their precise connectivity remains 
unclear. Thus, qualitatively at least, one could try to model the pacemaker 
cell network with the Kuramoto model. Now consider how this network might 
respond to an imposed cycle of light and dark (information of this sort is known 
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to be conveyed from the eyes to the pacemaker through a speciaHzed neural 
pathway). If the Ught-dark cycle is 24 hours long, we expect the electrical 
rhythms of many individual pacemaker cells to successfully entrain to it. But 
what if we alter the period or strength of the external forcing, as has been done 
in countless experiments on mice, hamsters, primates, and human volunteers 
[Moore-Ede et al. 1982]? Or what happens if the experiment is conducted 
on mutant organisms [Dunlap 1993, Takahashi 1995] whose intrinsic periods 
are a few hours longer or shorter than normal, or which may be intrinsically 
arrhythmic, having almost no free-running circadian rhythm at all? 

Questions like this can be addressed, in mathematically idealized form, 
within the framework of the periodically forced Kuramoto model [Sakaguchi 
1988, Antonsen et al. 2008, Ott and Antonsen 2008]. Its governing equations 
are given by 

^ = CO, + - J^ sm{'&j -^^)+F sm{at - i?,), (1) 

for i = 1, . . . ,N. Here t9j is the phase of oscillator i, cji is its natural frequency, 
K is the coupling strength, F is the forcing strength, a is the forcing frequency, 
and A^ ^ 1 is the number of oscillators. The natural frequencies are randomly 
distributed with a density g{uj), assumed unimodal and symmetric about its 
mean value ujq. 

This system is capable of rich dynamics because of its interplay among ran- 
domness, coupling, and forcing. The randomness comes from the variance of the 
natural frequencies. This effect tends to desynchronize the oscillators and scat- 
ter their phases. The coupling, on the other hand, tends to align the oscillators 
to the same phase, although it does not favor any particular frequency for the 
collective oscillation. In contrast, the forcing does favor a specific frequency, 
namely that of the external drive. Depending on the relative magnitudes of 



1 INTRODUCTION 

these competing effects, we expect to see various kinds of cooperative behavior 
and transitions between them. 

Before continuing, it proves useful to simphfy the governing equations in 
two ways. First, if we view the dynamics in a frame co-rotating with the drive, 
the exphcit time dependence in ([T|) disappears. To achieve this, let 

Oi = ^i- at. (2) 

Then ([U yields 

-^ = {oji - a) + -Y. sin(^j -ei)-F sin Oi, (3) 

for i = 1, . . . ,N. Second, as Kuramoto originally pointed out, it is helpful to 
introduce a complex order parameter z, given by 



Then the sum in ([3|) reduces to Im(i^ze^*^'), an identity which will prove useful 
later. 

The order parameter also has a nice physical interpretation. Its amplitude 
\z\ quantifies the phase coherence of the population: an incoherent state has 
z = 0; a perfectly coherent state has 1^1 = 1. Furthermore, the argument of z 
can be interpreted as the average phase of all the oscillators. So in a sense, the 
single complex number z{t) serves as a proxy for the state of the population as 
a whole. 

Sakaguchi [1988] was the first to study the periodically forced Kuramoto 
model. He derived a self-consistent equation for steady-state values of \z\, 
under the assumption that z{t) was entrained by the external force (meaning 
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that z{t) appeared motionless in the rotating frame). In numerical simulations 
of Eq. ([3]), however, Sakaguchi found that this state of "forced entrainment" 
was not always attained. For some values of the parameters, the system could 
settle instead into a state of "mutual entrainment." In this case a macroscopic 
fraction of the system self-synchronized at a different frequency from that of 
the drive, indicating that this sub-population had broken away and established 
its own collective rhythm. (For circadian rhythms, this would mean that the 
animal's internal clock was drifting relative to the outside world.) Sakaguchi's 
numerics further indicated how forced entrainment could be lost and give way 
to mutual entrainment. Such transitions were found to occur via two different 
mechanisms, corresponding to a pair of distinct bifurcation curves in parameter 
space. These curves appeared to join at a point, but Sakaguchi was unable to 
resolve the details of the cross-over region numerically. 

More recently, Antonsen et al. [2008] gave an improved analytical treat- 
ment of the model. Their linear stability analysis and numerical simulations 
also revealed an intriguing set of bifurcation curves, but the way the various 
curves join together still remained unclear. The overall layout of the stability 
diagram suggested that an underlying two-dimensional system was controlling 
the dynamics — a remarkable finding, given that the model ^ is essentially 
infinite-dimensional (recall A^ ^ 1). 

This tantalizing clue led Ott and Antonsen to an important discovery [Ott 
and Antonsen 2008]. They found that the Kuramoto model possesses an in- 
variant manifold, a special family of states for which the macroscopic dynamics 
becomes /ow-dimensional. In particular they showed that on this invariant man- 
ifold, the order parameter for the forced Kuramoto model ^ exactly satisfies 
a iwo-dimensional dynamical system, for the special case where the frequency 
distribution g{uj) is Lorentzian and the initial state satisfies certain strong an- 
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alyticity properties with respect to to. 

In this paper we analyze the two-dimensional system derived from the anal- 
ysis of Ott and Antonsen [2008]. Our results give the first complete picture 
of the bifurcation structure for the forced Kuramoto model. We obtain ex- 
plicit formulas for the system's saddle-node and Hopf bifurcation curves, as 
well as the codimension-2 Takens-Bogdanov point from which they emanate. 
Bifurcation theory predicts that a curve of homoclinic bifurcations should also 
emerge from the Takens-Bogdanov point; we compute this homoclinic curve 
numerically. 

The rest of the paper is organized as follows. Section 2 reviews the approach 
of Ott and Antonsen [2008] , leading up to their derivation of the reduced equa- 
tions for the order parameter dynamics. Section 3 presents new results about 
the bifurcations in this system and resolves the issue of how all the transition 
curves fit together. The final section discusses the implications of the results, 
their relation to prior work, the limitations of the approach used here, and some 
of the questions that remain. 

2 Derivation of the reduced equations 

The analysis of ([3]) is carried out in the continuum limit N ^ oo. Then the 
state of the system is described by a density function f{9, uj, t). Here / is defined 
such that at time t, the fraction of oscillators with phases between 9 and 9 + d9 
and natural frequencies between u; and uo + doj is given by f(9,uj, t) d9 duj. Thus 

/oo f2n 
/ f{9,oj,t)d9duj = l (5) 

-oo Jo 



and 

r2T7 

f{uj,9,t)d9 = g{uj), (6) 
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by definition of g{uj). 

The evolution of / is given by the continuity equation 



!->-■ 



(7) 



which expresses the conservation of oscillators of frequency lo. Here v{9, oj, t) is 
the velocity field on the circle corresponding to ^ as A^ ^ oo: 

/oo r2n 
/ sm{e' -6) f{e\uj',t)de'du' -F sine. (8) 
-oo Jo 

This expression can be written more compactly in terms of the complex order 
parameter z, which in the continuum limit becomes 



/oo /'27r 
/ e'^f{e,uj,t)dedu;. 
-co Jo 



(9) 



Using the identity mentioned in the Introduction, we note that the double 
integral in ([5]) simplifies to lTa{Kze~^^). Hence the continuity equation becomes 



dt dd^-' 



{iu-a) + -\{Kz + F)e-'^ - (Kz + F)*e^^| 



0, (10) 



where the asterisk denotes complex conjugation. 

Normally one would try to solve (jlOp by expanding / as a Fourier series in 



f{0,u;,t) 



2tt 



l + £A(a;,t)e*"^ + c.c. 



n=l 



(11) 



where c.c. denotes complex conjugate. Substitution of ([TT]) into ([9]) and ([TO] 
would generate an infinite set of coupled nonlinear ordinary differential equa- 
tions for the amplitudes fn{uj,t). Unfortunately the dynamics of this infinite- 
dimensional system would likely be difficult to analyze further. 
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It was at this point that Ott and Antonsen [2008] noticed something won- 
derful. They restricted attention to the special family of densities / for which 

fniu;,t) = [aiu;,t)r, (12) 

for all n > 1. In other words, they assumed that all the amplitudes fn are 
n}^ powers of the same function a{LO,t). Amazingly, this ansatz satisfies the 
amplitude equations for all n, so long as a evolves according to 

^ = ^{Kz + FT - i{u; - a)a - \{Kz + F)a' (13) 

and z satisfies 

POD 

z{t) = / a*{to,t)g{to)(ko. (14) 



Then, by further assuming that g{u) is a Lorentzian, 






and that a(w, t) satisfies certain analyticity conditions in the complex w-plane, 
Ott and Antonsen [2008] evaluated (J14p by contour integration and thereby 
derived the following exact evolution equation for the order parameter z: 



d? 1 

- = - [{Kz + F)- {Kz + FTz^] - [A + i{a - u;o)] z. (16) 



The conditions required were that a{uj,t) can be analytically continued from 
the real w-axis into the lower half of the complex w-plane for all t > 0; that 
|a(a;,t)| ^ as Im(a;) — > — oo; and that |a(a;,0)| < 1 for real io. 



3 ANALYSIS OF THE REDUCED EQUATIONS 11 

3 Analysis of the reduced equations 
3.1 Scaling the equations 

We turn now to the analysis of the two-dimensional system (|16p . The first step 
is to reduce the number of parameters by nondimensionalizing the system. Let 
i = At,F = F/A, k = K/A, a = a/ A and ujq = ujq/A. Then the form of ^ 
stays the same except that A no longer appears (in effect, A has been set to 
1 without loss of generality) and all the other parameters now have hats over 
them. For ease of notation we drop the hats in what follows, remembering that 
all the parameters are now dimensionless. Also, let 

VL = a — ojQ (17) 

denote the (dimensionless) detuning between the drive frequency and the pop- 
ulation's mean natural frequency. Then if we introduce polar coordinates 

z = pe''^ (18) 

and separate (J16p into real and imaginary parts, we obtain the dimensionless 
evolution equations for p and (j): 

p' = ^p(l-p^)-p + ^{l-p^)cos4> (19) 

(p' = - n + ^(p+-jsmcf) (20) 

where the prime denotes differentiation with respect to dimensionless time. 
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3.2 Stability diagram and phase portraits 

Our next goal is to obtain the stability diagram for Eqs. P^ - ([2U|) . Before 
delving into the details, which can become intricate at times, we jump to the 
final result. Figure 1 shows the stability diagram for the representative case 
where K = 5. Here the various stability regions labeled A-E correspond to the 
phase portraits shown in Fig. 2. 

We realize that these figures appear complicated at first glance, so let us be- 
gin by offering a few general remarks about them. Figure 1 is divided into five 
regions, A-E, by the bifurcation curves labeled saddle-node, Hopf, homoclinic, 
and SNIPER. In the places where two or more of these curves nearly coincide, 
Fig. 1(a) becomes especially confusing. To clarify what is going on in such re- 
gions, Figs. 1(b) and 1(c) zoom in near two codimension-2 points of interest (to 
be discussed in detail later). Since even these figures can be hard to interpret, 
we have tried to make everything as clear as possible by presenting a schematic 
Fig. 1(d). Unlike Figs. l(a)-(c), which are numerically accurate. Fig. 1(d) is 
only topologically correct. We have distorted some of stability regions and 
pulled certain curves apart to make the layout of the diagram transparent, and 
to highlight the three different codimension-2 points that will later be seen to 
organize the entire diagram. 

A similar but incomplete version of Fig. 1 was obtained previously by An- 
tonsen et al. [2008]; see Fig. 3 in their paper. Those authors generated their 
results based on direct simulations of Eq. Q for N = 1000 oscillators. They 
also compared their numerics to analytical results they derived for the existence 
and stability of equilibrium points for ([3]), which correspond to entrained states 
in the original frame. Our approach, in contrast, is to analyze the reduced 
system Eqs. ()19p -()20p. We do not present numerical results for the full system 
([3]) because in every case we have checked, our results match those reported 
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already by Antonsen et al. [2008], except in cases where the previous methods 
were inconclusive. 

3.3 Saddle-node and SNIPER bifurcations 

It is algebraically awkward to solve for the fixed points of Eqs. (I19p - (j20p in 
terms of the parameters. Fortunately, we do not need to solve for them. Since 
we are mainly interested in the bifurcation curves, we can make headway more 
easily by imposing an appropriate bifurcation condition and then solving for 
the parameters in terms of the fixed point, rather than the other way around. 
This is a standard trick in bifurcation theory, and it allows us to derive the 
bifurcation curves in closed form, either explicitly or as parametric equations. 

For example, at a saddle-node bifurcation, one of the eigenvalues equals 
and hence the determinant of the Jacobian vanishes there. (The same would 
be true at transcritical or pitchfork bifurcations, but given the absence of the 
constraints or symmetries associated with these types of bifurcations, there's 
no reason to expect either of them to occur here.) 

Hence to find the locus of saddle-node bifurcations, we solve /?' = 0, </>' = 
and 6 = simultaneously, where S denotes the determinant of the Jacobian. 
The trick is to regard the unknown values of the variables p and (p on equal 
footing with the parameters K, $7 and F. Then the resulting system of 3 equa- 
tions in 5 unknowns can be solved explicitly to yield a parametrization of the 
saddle-node bifurcation surface. Various parametrizations are possible. One 
convenient choice is to express the parameters in terms of the fixed-point val- 
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ues of p and </>. We find that the saddle-node surface is then given by 
^ ^ 2{p^ + 2p^cos2<p+l) 

(l-p2)2(l + p2cos20) 

^ ^ {p^ + p)^ sm2(t> 

(l-p2)2(l + p2cos20) 



4/9^ (p2 + l) cos (j) 
(l-p2)2(l+p2(,os20) 



(23) 



where we allow p and </> to sweep over their full ranges < p < 1, — vr < (j) < tt. 

This parametrization provides some interesting information. For instance, 
it shows that K increases monotonically with p, for each fixed value of (j). Hence 
K >2 everywhere on the saddle-node surface, with the minimum value K = 2 
being attained when /? = and hence F = 0, or in other words, when there 
is no forcing. This result makes sense. In the absence of forcing, the system 
is just the traditional Kuramoto model with a Lorentzian g{io), and K = 2A 
(or in dimensionless terms, K = 2) \s the well-known formula for the critical 
coupling at the onset of mutual synchronization [Kuramoto 1975, 1984]. 

To compare our results with those obtained numerically by Antonsen et al. 
[2008] , it is more illuminating to slice through the saddle-node surface at a fixed 
value of i<' > 2 and then plot the resulting saddle-node curves in the (fi, F) 
plane. To find these curves we solve p' = 0, (/>' = and 5 = for i7,sin(/> and 
cos (j), and then use sin (f) + cos2 = 1 to solve for F, now regarding K and p 
as parameters. The result is the following parametrization of the saddle-node 
curve: 



[p^ + lf\lK{p^ -1){k{p^ -if -A] -A{p^ + 1) 
2(/>2_i) 



^SN = '- T-v— \2 (24) 



y/2p^jK'^ (p2 - 1)3 + 2K (p4 _ 4p2 + 3) _ 8 
i^SN = ^^^—^, . (25) 
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Figure 1 plots this saddle-node curve for the case K = 5, as previously studied 
by Antonsen et al. [2008]. We compute the curve for all values < p < 1, 
disregarding any values that yield non-real results for fi or F. 

The two branches of the saddle-node curve intersect tangentially at a codimension- 
2 cusp point, as highlighted in Fig. 1(c) and marked schematically in Fig. 1(d) 
by the solid square . For K = 5, the parameter values at the cusp are ~ 3.5445 
and F « 3.4164. 

Along with local saddle-node bifurcations, the lower branch of the saddle- 
node curve(where F ~ J7) also includes a large section consisting of saddle- 
node infinite-period (SNIPER) bifurcations. These have important global im- 
plications, because they create or destroy limit cycles in the phase portrait of 
Eqs. (HSD-dSni). 

3.4 Hopf bifurcation 

Next we calculate the locus of parameter values at which Hopf bifurcations 
occur. We impose the fixed point conditions (j)' = 0, p' = as before, but 
now require that the Jacobian has zero trace and positive determinant — the 
latter two conditions are equivalent to requiring that the eigenvalues be pure 
imaginary. 

Solving simultaneously for (p' = 0, p' = and trace = 0, we find 

sine/) = , ^ — ' (26) 

F^/I{^^K^/I{T2 

COS,/. = _ 1:^:1^ (27) 

2F^/I{T2 

Because p depends only on K, we can go a bit further than we did in the 
saddle-node case. Using sin^ -|- cos^ </> = 1 as before, F can now be obtained 
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explicitly in terms of K and 0: 



i^Hopf-^y irr2 ^^^^ 

For the special case K = 5 studied by Antonsen et al. [2008], Eq. ([29|) becomes 

i^Hopf = 3^ Y 7 ^225 + 19602 (30) 

Figure 1 plots the graph of -Ffjopfl^)- Notice how straight it is, and that it 
nearly lines up with the lower branch of the saddle-node curve. 

3.5 Takens-Bogdanov point 

As mentioned above, for Eq. (j29p to truly signify a Hopf bifurcation, the Ja- 
cobian determinant must be positive at the corresponding parameter values 
(fi, F) . This will be the case if 0, and F are sufficiently large. Specifically, 
their values must exceed those at the Takens-Bogdanov point 



„,, . iJi^^ (31) 



4:{K + 2) 


\k 2)J^'- 


-2K^ + AK -8 


4^^ ^^V 


K + 2 



Ftb = jiK-2)^ ^^^ (32) 

obtained by solving four simultaneous equations: (p' = 0, p' = 0, trace = 0, and 
determinant = 0. 

The Takens-Boganov point is marked with a filled circle on Figs. 1(a) and 
1(d). In addition to serving as the endpoint of the Hopf curve, it splits the 
upper branch of the saddle-node curve into two sections of different dynamical 
character. On the lower section (below the Takens-Boganov point), an unsta- 
ble node collides with a saddle along the saddle-node bifurcation curve; this 
can be seen by comparing regions D and A, as shown in Figs. 2(d) and 2(a). 
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The opposite situation occurs on the upper section of the saddle-node curve, 
above the Takens-Boganov point. Here a stable node collides with a saddle, 
corresponding to the transition between regions B and A; see Figs. 2(b) and 
2(a). 

3.6 Homoclinic bifurcation 

The theory of the Takens-Bogdanov bifurcation implies that a curve of homo- 
clinic bifurcations must also emerge from the codimension-2 point, tangential 
to the saddle-node and Hopf curves. For the case K = 5 shown in Fig. 1, 
JItb = 28 ^^'^ -^TB = j\/^- The curve shown in the diagram was com- 
puted numerically. It is almost indistinguishable from the Hopf curve and thus 
produces a very small region between them, as shown in Fig. 1(b). 

A striking feature of the homoclinic curve is that after moving parallel to 
the Hopf curve for a while, it makes a sharp backward turn and then joins 
onto the lower branch of the saddle-node/SNIPER curve, meeting that curve 
tangentially at a codimension-2 "saddle-node-loop" point [Guckenheimer 1986, 
Izhikevich 2000] marked by a filled diamond in Figs. 1(b) and 1(d). 

3.7 Phase portraits and bifurcation scenarios 

As we have seen, the bifurcation curves in Fig. 1 partition the stability diagram 
into five regions, labeled A-E. We now give a fuller treatment of the dynamics 
associated with each region and the transitions from one to another. 

3.7.1 Region A: Forced entrainment 

Here the order parameter z approaches a stable fixed point for all initial condi- 
tions, as shown in Fig. 2(a). To interpret what this means, recall that all our 
analysis has assumed a frame co-rotating with the drive. Hence this stable fixed 
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point represents a state in which the order parameter is moving periodically 
while staying phase-locked to the drive. Therefore, back in the original frame, 
a macroscopic fraction of the oscillator population must also be moving in rigid 
synchrony, locked to the same frequency as the drive signal. 

3.7.2 Region B: Bistability between two states of forced en- 
trainment 

Now suppose we weaken the forcing. Imagine moving down along a vertical 
line in Fig. 1(b), decreasing F while holding fi fixed. As we do so, we first 
pass from region A into the extremely narrow region B by crossing through the 
upper branch of the saddle- node curve (j24p . At this bifurcation, a stable node 
is born out of the vacuum, along with a saddle point. Meanwhile, the stable 
fixed point that we encountered in Region A still exists; it lies in the lower right 
part of Fig. 2(b). 

Thus Region B depicts a form of bistability. Depending on the initial con- 
ditions, the system chooses one of two possible states of forced entrainment, 
differing in their phase coherence (the magnitude of z) and their phase rela- 
tionship to the drive signal (the argument of z). 

3.7.3 Region C: Bistability between forced entrainment and 
phase trapping 

Continuing our vertical descent through Fig. 1(b), we next cross from B into C 
by passing through the curve of Hopf bifurcations, Eq. (j29p . The stable fixed 
point created in Region B now loses stability and gives birth to a tiny attracting 
limit cycle (Fig. 2(c)). On this cycle the order parameter still runs at the same 
average frequency as the drive but its relative phase and amplitude now wobble 
slightly. Because these variations remain trapped within tight limits, one says 
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the system is phase trapped (as opposed to phase locked) to the drive. Back in 
the original non-rotating frame, the macroscopic dynamics for this state would 
be quasiperiodic with two frequencies. This is not the only attractor, of course; 
the state of forced entrainment seen earlier in A and B persists, so we still have 
bistability, but now between phase trapping and forced entrainment. 

3.7.4 Region D: Forced entrainment 

Passing from Region C to D carries us across a curve of homoclinic bifurcations. 
As we approach this curve from above, the tiny limit cycle in Fig. 2(c) expands. 
At the bifurcation it touches the saddle point and forms a homoclinic orbit. 
Beyond the bifurcation the phase portrait looks like that shown in Fig. 2(d). 
An invariant loop has been created, in which the saddle and the original stable 
node are now connected by both branches of the saddle's unstable manifold. 
The stable node is the unique attractor. Hence the system again falls into a 
state of forced entrainment. 

3.7.5 Region E: Mutual entrainment 

Forced entrainment is finally lost when we pass from Region D to E. When 
crossing the lower branch of the saddle-node curve, we need to be careful to 
specify exactly where the crossing occurs. Specifically, do we cross to the left or 
right of the codimension-2 saddle-node-loop point (filled diamond in Fig. 1(b)) 
at which the homoclinic curve meets the saddle-node curve? 

Suppose first that we cross below and to the left of the saddle-node-loop 
point. Then in Fig. 2(d) the saddle and node would slide toward each other 
along the invariant loop, coalesce, and disappear, leaving a stable limit cycle in 
their wake. Thus, this saddle-node bifurcation is actually a SNIPER (saddle- 
node infinite-period) bifurcation. 
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The limit cycle created by the bifurcation is globally attracting. Hence the 
order parameter always settles into periodic motion in the rotating frame. But 
unlike the limit cycle of Fig. 2(b) this cycle winds around the origin of the 
z-plane, marked by an asterisk in Fig. 2(d). This is an important distinction, 
because it implies that the phase of z now increases monotonically relative to 
that of the drive. Consequently the order parameter z{t) oscillates at a different 
average frequency from the drive signal, implying that a macroscopic fraction 
of the oscillator population has broken loose from the drive. In other words, 
the system has spontaneously mutually entrained itself, at least in part. This is 
therefore one mechanism by which forced entrainment can give way to mutual 
entrainment. 

But there are other possible mechanisms as well. For example, consider 
Fig. 1(b) again, and now direct your attention to the upper right corner. By 
moving down along the right side of the picture, we can cross directly from 
C to E, without ever going through D. This happens when we cross through 
the portion of the lower saddle-node curve lying above and to the right of the 
saddle-node-loop point. In this case the bifurcation is not a SNIPER; it's just 
an ordinary saddle-node bifurcation. To visualize this scenario, imagine sliding 
the saddle in the middle of Fig. 2(c) to the right along its unstable manifold 
until it collides with the node and annihilates it. During this process the limit 
cycle in Fig. 2(c) grows. And so the phase portrait now resembles the one 
shown in Fig. 2(e). 

A third scenario is much simpler. Suppose fi > J^cusp) so that we're well to 
the right of the cusp in Figs. 1(c) and 1(d). Then as we decrease F, we move 
directly from A to E. Forced entrainment gives way to mutual entrainment 
through a supercritical Hopf bifurcation. 
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4 Discussion 

4.1 Stability diagram 

The main result of the paper is the stabihty diagram shown in Fig. 1. We have 
focused on the analytical derivation of several of the curves in this picture and 
tried to clarify how they fit together and what they imply about the system's 
overall dynamics. Having immersed ourselves in the details, it is worthwhile to 
step back and try to understand the broader lessons that this picture holds. 

Figure 1 essentially divides into two big regions. One represents forced 
entrainment, wherein a macroscopic fraction of the population is phase-locked 
to the drive. The rest of the population consists of oscillators belonging to 
the infinite tails of the frequency distribution; these remain unlocked. Thus it 
would be more accurate to speak of "partial forced entrainment," though we 
hope the intended meaning of the shorter name is clear. 

The other main region represents (partial) mutual entrainment. Now there 
are three qualitatively different groups of oscillators: (1) the unlocked oscillators 
in the tails; (2) the oscillators entrained by the forcing; and (3) a self-organizing 
group of oscillators that entrain one another at a frequency different from that 
of the drive. The existence of this third group causes the order parameter to 
wobble or drift periodically relative to the drive, as manifested by a stable limit 
cycle in the phase portraits (Figs. 2(c) and 2(e)). 

4.2 Comparison to Adler equation 

The boundary between forced and mutual entrainment is complicated when 
viewed at a fine scale, as shown in Fig. 1(b). But from a bird's-eye view, it 
looks very much like the straight line F = 0. Here's why: this is the result one 
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would expect from the Adler equation 

(j)' = -n-Fsm(l>, (33) 

which has been used to model the entrainment dynamics of phase-locked loops 
[Adler 1946], lasers [Siegman 1986, Yeung and Strogatz 1998], and fireflies [Er- 
mentrout and Rinzel 1984], among many other systems. The two-dimensional 
system (fT9]) - (pO]) reduces to Adler's equation as ET ^ oo, in the sense that p 
approaches 1 on a fast time scale, while (p obeys (|33p on a slow time scale. 

The intuitive explanation is that in this limit, the coupling between oscilla- 
tors is so strong that the population acts like one giant oscillator, with nearly 
all the microscopic oscillators at the same phase. Hence the order parameter 
amplitude remains close to /> = 1 at all times, so the system behaves as if it 
had a very strongly attracting limit cycle. This explains why the dynamics of 
the forced Kuramoto model mimic the Adler equation in this limit. 

For a more analytical route to the same conclusion, look at the large--ftr 
behavior of the Takens-Bogdanov point, which essentially lies on the dividing 
line behind the two big regions. The formulas (J3ip - (j32p imply that 

^ ~ 1 - 8if-^ (34) 

as i^ — > OO. Thus F ~ il for large and even moderate values of K. 

4.3 Comparison to forced van der Pol equation 

For weaker coupling, but still large enough that the system can partially self- 
synchronize {2 < K < oo), the population again acts like a single limit cycle 
oscillator, but now with a limit cycle that is only weakly attracting. As before, 
the complex order parameter plays the role of this effective limit-cycle oscillator. 
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So when forcing is applied, we expect the overall dynamics to be like those 
of a forced, weakly nonlinear oscillator. And indeed, the stability diagram 
bears a striking resemblance to that of a forced van der Pol oscillator, in the 
limit of weak nonlinearity, weak detuning, and weak forcing. As in the problem 
studied here, the stability diagram for this well-studied system [Guckenheimer 
and Holmes 1983] is also organized around a Takens-Bogdanov bifurcation and 
a saddle-node-loop bifurcation. 

Likewise, some of the regions in the van der Pol diagram are unusually thin 
and small. This helps to explain why they were overlooked for decades, until 
the theory of the Takens-Bogdanov bifurcation was developed and guided later 
researchers to the missing transitions that, on topological grounds, had to be 
there. 

One always expects small regions in systems with Takens-Bogdanov bifur- 
cations because, according to normal form theory, the saddle-node, Hopf, and 
homoclinic curves have to intersect tangentially at the Takens-Bogdanov point. 
But here, as in the van der Pol problem, the regions are even smaller still, 
because they must also hug the line F ~ $7, for the reasons given above. 

4.4 Caveats 

It is important to understand what has — and has not — been shown by the anal- 
ysis presented in this paper. Following Ott and Antonsen [2008], we made 
a number of very particular choices in the course of reducing an infinite- 
dimensional problem to a two-dimensional one. We chose a special family of 
initial states (see Eq. p2|) ) and showed that they formed an invariant manifold. 
In other words, if the condition (J12p is satisfied initially, it is automatically satis- 
fied for all time. Then we chose a special distribution of natural frequencies (see 
Eq. ()15p ). and required further that the initial state satisfies certain strong an- 
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alyticity properties with respect to its dependence on these frequencies. Taken 
together, these choices then imphed that the system's order parameter evolves 
according to the two-dimensional dynamical system (|16p . 

If the conclusions that followed were sensitive to these choices, we would 
not have accomplished much. But there is reason to believe that the results 
are robust, and largely independent of these choices. The strongest evidence is 
numerical. Every time we have run simulations of the forced Kuramoto model 
([1]) for hundreds or thousands of oscillators, we have seen all the attractors and 
bifurcations predicted by the analysis, where they are supposed to be. Ott and 
Antonsen [2008] found similar agreement when they studied other variants of 
the Kuramoto model. 

This suggests that the flow on the invariant manifold faithfully captures 
the macroscopic dynamics of the full system, at least in some sense. Unfor- 
tunately, we do not know how to make this statement precise. The issue is 
probably subtle. We do not believe, for example, that the invariant manifold 
is everywhere transversely attracting — it certainly isn't in other problems we 
have studied. For example, applying the method of Ott and Antonsen [2008] 
to the Kuramoto model with a bimodal frequency distribution, we found that 
the invariant manifold in that case could be transversely repelling at certain 
points [Martens et al. 2008]. 

Nor are we sure whether all the attractors for the full system lie within the 
invariant manifold. If they did, that would explain why this manifold controls 
the system's long-term macroscopic dynamics. But we have no proof of this 
weaker statement either. 

Now regarding the choice of a Lorentzian frequency distribution: this was 
crucial to the analysis, but not, we suspect, to the results. Sakaguchi [1988] 
used a Gaussian g{uj) and found the same attractors and bifurcations as we 



5 ACKNOWLEDGEMENTS 25 

did. Our own simulations for the Gaussian case (unpublished) show that the 
stability diagram is different in numerical details, of course, but its topology is 
unaffected. 

On the other hand, the algebraic form of the forced Kuramoto model, with 
its purely sinusoidal coupling and driving, probably is crucial. The ansatz (|12|) 
no longer works if the model contains higher harmonics. Indeed, the bifurcation 
behavior of the classical (unforced) Kuramoto model is known to be altered 
when generic periodic functions are used in place of a pure sine function in the 
coupling [Daido 1994, Crawford 1995]. So we expect new phenomena to arise 
in the forced Kuramoto model as well, when one departs from pure sinusoids 
in the driving and coupling terms. 
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Figure Captions 

Figure 1: Stability diagram for the forced Kuramoto model. Bifurcation 
curves are shown with respect to the strength F and detuning fi of the external 
forcing, both of which have been non-dimensionalized by the width A of the 
distribution of the oscillators' natural frequencies. The dimensionless coupling 
strength is fixed at K = 5. 

(a) Regions A-E correspond to qualitatively different phase portraits; see 
text and Fig. 2 for explanations. Four types of bifurcations occur: supercritical 
Hopf bifurcation; homoclinic bifurcation; and two types of saddle-node bifurca- 
tions. The saddle-node bifurcations on the upper branch, and also those on the 
lower branch between the cusp and the saddle- node- loop point, are purely local. 
In contrast, those on the portion of the lower branch extending from the origin 
to the saddle-node-loop point have global consequences; they are saddle-node 
infinite-period bifurcations, or SNIPERs, which create or destroy limit cycles. 
The filled circle marks a codimension-2 Takens-Bogdanov point, at which the 
Hopf, homoclinic, and upper saddle-node curve intersect tangentially. 

(b) Enlargement of the cross-over region, just to the right of the Takens- 
Bogdanov point, where all four bifurcation curves run nearly parallel to one 
another. Three of them (Hopf, SNIPER, and the lower branch of saddle-node 
bifurcations) meet at a codimension-2 saddle-node-loop point, marked by a 
filled diamond. 

(c) Enlargement of the region near the codimension-2 cusp point (filled 
square), where the upper and lower branches of saddle- node bifurcations meet 
tangentially. The two branches are almost indistinguishable in this image. 

(d) Schematic version of the stability diagram, intended to show how the 
bifurcation curves connect in the confusing cross-over region. Tangential inter- 
sections have been opened up for clarity. 
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Figure 2: Phase portraits for the two-dimensional dynamics of the com- 
plex order parameter z, or equivalently, for the variables p, (p regarded as polar 
coordinates. Open dots, unstable fixed points. Closed dots, stable fixed points. 
Asterisk, origin of the z-plane. Dashed curves, stable and unstable manifolds 
of the saddle point. The panels are not all shown at the same scale; the regions 
shown in (b) and (c) are small and have been blown up here for clarity. 



